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ABSTRACT 

We demonstrate that the nuclei of galaxies containing supermassive black holes can be triaxial in shape. 
Schwarzschild's method was first used to construct self-consistent orbital superpositions representing 
nuclei with axis ratios of 1 : 0.79 : 0.5 and containing a central point mass representing a black hole. 
Two different density laws were considered, p cx r -7 , 7 = {1, 2}. We constructed two solutions for each 7: 
one containing only regular orbits and the other containing both regular and chaotic orbits. Monte-Carlo 
realizations of the models were then advanced in time using an Y-body code to verify their stability. 
All four models were found to retain their triaxial shapes for many crossing times. The possibility that 
galactic nuclei may be triaxial complicates the interpretation of stellar-kinematical data from the centers 
of galaxies and may alter the inferred interaction rates between stars and supermassive black holes. 
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1. INTRODUCTION 

The phenomenon of triaxiality has remained of central 
importance to our understanding of galaxy dynamics since 
the demonstration by Schwarzschild (1979) of the exis- 
tence of self-consistent triaxial equilibria. Schwarzschild's 
models had large, constant-density cores and the major- 
ity of the orbits were regular, respecting three integrals of 
the motion. But the realization that the central densities 



of elliptic al galaxies and bulges are very high ( Crane et 
al. 1993), and that supe rmassive black holes are generi c 
components of galaxies ( Kormcndy fc Richstonc 1995| ), 



has modified somewhat our ideas about the persistence 
of triaxiality. A central black hole can destroy triaxiality 
on large sc ales by rendering the cen ter-filling box orbits 
stochastic ( Gerhard fc Binney 1985 ). Evolution to glob- 
ally axisymmetric shapes can occur in just a few cross- 
ing times if th e black hole contains of order 10~ 2 of t he 
galaxy's mass ( [Merritt fc Quinlan 19981 ; [Bcllwood 200l] ). 

Less is known about the possibility of maintaining triax- 
ial shapes at the very centers of galaxies, where the grav- 
itational force is contributed roughly equally by the stars 
and by the nuclear black hole. Regular orbits, both box- 
like and tube-like, exist in this region (iMerritt fc Vallurj 
1999[ |Sambhus fc Sridhar 2000| ; [Poem fc Merritt 2001| , Pa- 
per I) but the box- like orbits (the "pyramids" ) have shapes 
that are generally contrary to the figure elongation, mak- 
ing them less useful than classical box orbits for main- 
taining a triaxial shape. Furthermore the pyramid orbits 
disappear within the "zone of chaos," which begins at a 
radius where the enclosed stellar mass is a few times that 
of the black hole (Paper I), roughly 10 2 pc in a typical 
bright elliptical galaxy. 

Here we demonstrate that triaxial equilibria do in fact 
exist in the vicinity of the black hole. We first (§2) use 
Schwarzschild's technique to construct self-consistent su- 
perpositions of orbits computed in a fixed triaxial poten- 
tial representing the stars and a central point mass. We 
then (§3) test the long-term stability of the models by us- 
ing an Y-body code to advance them forward in time; we 



find that the Y-body models maintain their triaxial shapes 
for many crossing times. Finally we present some of the 
observable properties of the models (§4) and discuss the 
implications of triaxiality for observational and theoretical 
studies of galactic nuclei (§5). 

2. SCHWARZSCHILD SOLUTIONS 

We model the stellar nucleus as a triaxial spheroid with 
a power-law dependence of density on radius, 

(la) 

I (lb) 



p± = p m 



y 



inside of some surface m — m ou t discussed below. We took 
7 = 1 or 2; these values correspond roughly to the density 
profiles observed at the ce nters of bright and fai nt ellipti- 
cal galaxies, respectively ( |Gcbhardt et al. 1996|) . The 
triaxiality T is defined as 



T 



(2) 



We chose a = 1.0,6 = 0.79, c = 0.5 for both models, cor- 
responding to T — 0.5, "maximal" triaxiality. The central 
black hole is represented by a point mass with M^ = 1. 
Expressions for the gravitational potential and forces cor- 
responding to this mass model are given in Paper I. 

Although our mass model for the stars is scale-free, the 
presence of the black hole imposes a scale. We identify 
three characteristic radii associated with the black hole. 
At a radius of r g , the enclosed mass in stars (defined as 
the mass within an equidensity surface which intersects 
the x-axis at x — r g ) is equal to that of the black hole. 
In real galaxies, r g is typically a factor of ~ a few greater 
than the "sphere of influence" 77,., where 



GM, 



hi, 



10.8 pc 



/ Mbh 
\10 S M Q ) V200 km s 



(3) 



depends on the stellar velocity dispersion cr» The third 
fiducial radius is r c /,, the radius at which the character of 



2 



the box-like orbits undergoes a sudden transition to chaos 
in triaxial potentials. Table 1 gives values of r g and r c h for 
our two mass models; values of r c h are approximate and 
are taken from Paper I. Table 1 also gives values of the 
dynamical times at r g and r c h , defined as the period of a 
circular orbit of the same energy in the equivalent spherical 
potential, defined to have a scale length a ave = (abc) 1 ^ 3 . 
The enclosed stellar mass within r c h is about 3M;,/,, for 
7=1 and 6M^ for 7 = 2. 

Table 1 



7 = 2 





7=1 


7 = 2 


r 9 


0.64 


0.20 


T D (r g ) 


0.98 


0.17 


Teh 


1.1 


1.3 


T D (r ch ) 


1.8 


1.5 


m out 


5.6 


3.8 


TD(m out ) 


5.1 


4.6 



Any numerical realization of a stellar system must have 
an outer boundary. Our goal was to test self-consistency 
out to a radius of at least r c h, corresponding to ~ 2r g 
for 7 = 1 and ~ 6r g for 7 = 2. We therefore chose 
the outer surface of our mass model, m — m out , to be 
large enough that almost all of the density at r c h in a 
real galaxy would be contributed by orbits with apocen- 
ters below this surface. In order to estimate m ou t, the 
isotropic distribution functions corresponding to a spheri- 
cal galaxy with the density law of equation (1), 7 = {1,2}, 
and a central point mass were computed and transformed 
to F(r apo ;r), the distribution of apocentric radii r apo of 
orbits passing through r. For 7 = 1, we found that orbits 



with r apo < 5r c h contribute ~ 70% of the density at r c h- 
For 7 = 2, setting r apo w 3r c h is sufficient to give ~ 75% 
of the density at r c h- {r c h in the spherical models was 
defined as the radius containing the same mass as in the 
triaxial models.) Our corresponding choices for m ou t are 
given in Table 1. 

We followed standard procedures for cons tr ucting the 
Schw arzschild solutions flSchwarzschild 1993 ; Mcrritt &| 
Fridman 1996). The mass model within m out was divided 
into 64 equidensity shells and each shell was divided into 
48 cells per octant, for a total of 3072 cells. Shells were 
more closely spaced near the center (Figure 1). Orbits 
were computed in two initial condition spaces: stationary 
start space, which yields box-like orbits, and X — Z start 
space, which yields mostly tube orbits. Orbital energies 
were selected from a grid of 42(52) values for 7 = 1(2), 
defined as the energies of equipotential surfaces that were 
spaced similarly in radius to the equidensity shells. The 
outermost energy shell intersected the x axis at x — m ou t 
for both models. Orbits were integrated for 100 dynami- 
cal times, as defined above, and their contributions to the 
masses in the cells were recorded. In order to distinguish 
regular from stochastic trajectories, the largest Liapunov 
exponent was computed for each orbit. For 7 = 1, the 
total number of orbits integrated was 18144 of which 9751 
were found to be regular. For 7 = 2, the numbers were 
22464 and 15048 respectively. Orbital weights that repro- 
duced the masses in the 3072 cells w ere found via quadratic 
programming (e.g. Dcjonghc 1989| ). 




0.55 0.78 0.87 



Fig. 1. — Cumulative mass fractions F contributed by different 
kinds of orbits to different shells of the triaxial models. Box-like 
orbits are green, tube orbits (both 2- and ic-tubes) are orange, and 
chaotic orbits are purple. Higher energies are indicated via darker 
shades; the color bar relates the shade to the radius of the equipo- 
tential surface, as explained in the text. Numbers below the color 
bar indicate radii where the equidensity shells intersect the i-axis. 
(a), (c) correspond to solutions with only regular orbits for 7 = 1,2 
respectively, (b), (d) correspond to solutions with both regular and 
chaotic orbits. 

As expected, orbital superpositions that exactly repro- 
duced the cell masses in all of the cells, including the out- 
ermost ones, could not be found since only a few orbits 
visit the outer shells. We were able to find "exact" solu- 
tions (solutions which precisely reproduced the densities in 
a subset of cells) by relaxing the constraints in a number of 
ways; for instance, by eliminating the angular constraints 
on the outermost shells (forcing only the integrated shell 
mass to be fit), or by ignoring the outer cells entirely when 
imposing the constraints. For example, when only regular 
orbits were included, "exact" solutions could be found for 
7=1 for all cells within shell 51, corresponding to a ra- 
dius of ~ 3.2. This radius is substantially greater than r c h 
(Table 1). 

When these "exact" solutions were advanced forward 
in time, as described below, they were found to exhibit 
significant evolution at large radii, due to the fact that 
the model's density was not correctly reproduced in the 
outermost cells. We therefore constructed new solutions 
which were not "exact," but which were constrained to re- 
produce the densities in all cells within m ou t with as high 
an accuracy as possible. Such solutions exhibited smaller 
fractional errors in the innermost cells (r < r c h) than in 
the outer ones. Models constructed in this way were found 
to evolve much less than the "exact" solutions and provide 
the basis for the discussion below. 
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Table 2 



7=1 



7 = 2 





z-tubcs 


ir-tubes 


pyramids 


chaotic 


7=1: (regular) 










r < 0.5 


0.61 


0.10 


0.29 


0.00 


r < 1.0 


0.56 


0.08 


0.37 


0.00 


r < 1.5 


0.53 


0.07 


0.40 


0.00 


7 = 1: (all) 










r < 0.5 


0.09 


0.04 


0.26 


0.61 


r < 1.0 


0.10 


0.03 


0.25 


0.62 


r < 1.5 


0.10 


0.02 


0.25 


0.63 


7 = 2: (regular) 










r < 0.4 


0.75 


0.10 


0.15 


0.00 


r < 0.8 


0.73 


0.13 


0.14 


0.00 


r < 1.2 


0.73 


0.13 


0.13 


0.00 


7 = 2: (all) 










r < 0.4 


0.06 


0.05 


0.45 


0.45 


r < 0.8 


0.05 


0.05 


0.44 


0.47 


r < 1.2 


0.04 


0.05 


0.42 


0.49 



The orbital content of the latter solutions is described 
in Figure 1 and Table 2. We constructed two types of solu- 
tion for each 7: solutions in which only regular orbits were 
used; and solutions in which all orbits (regular and chaotic) 
were provided to the Schwarzschild algorithm. The domi- 
nant family of orbits in the models containing only regular 
trajectories was found to be the z-tubes, orbits which cir- 
culate about the short axis of the figure. Plots of the 
individual orbits used in the solutions confirm that many 
of the z-tubes have the correct shape for reproducing the 
triaxial figure, being more elongated in x (the long axis) 
than in y (the intermediate axis). Most of the remain- 
ing contribution to the density of the regular nuclei was 
found to come from pyramid orbits and high-energy box 
orbits in the 7=1 model, while for 7 = 2, roughly equal 
contributions came from pyramid orbits and from the x- 
tubes. (See Paper I, Figure 1 for illustrations of the var- 
ious orbit families.) The most heavily occupied pyramid 
orbits had shapes similar to "saucer" orbits, z-tubes asso- 
ciated with a 2 : 1 resonance (Paper I, Figure lc); both 
models contained roughly equal numbers of saucers and 
pyramids. High-energy box-like orbits are more impor- 
tant in the 7 = 1 case than in the 7 = 2 case, because: (1) 
high-energy tube orbits in the 7=1 model are elongated 
opposite to the figure; (2) most of the high-energy regular 
box orbits in the 7 = 2 model are bananas, which are not 
very useful for filling the outermost shells. 

In the models constructed using both regular and 
chaotic trajectories, the number of stars on chaotic orbits 
is surprisingly high: ~ 60% for 7 = 1 and ~ 45% for 7 = 2 
at all radii. Plots of the most highly populated chaotic or- 
bits suggest that they have reached a nearly steady state 
in a time-averaged sense, filling the volume defined by the 
equipotential surface corresponding to their energy. Sta- 
tionary chaotic building blocks like these have been used 
befor e when c onstructing triaxial models (e.g. [Merritt & 



Fridman 1996) but no published triaxial models have con- 
tained nearly so large fraction of stars on chaotic orbits. 
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Fig. 2. — Evolution of the ./V-body axis ratios; a = 1. r defines 
the longest dimension of the ellipsoid within which the axis ratios 
were determined, using the algorithm described in the text. Tpj is 
the dynamical time at this radius. Black (red) curves correspond to 
solutions where only regular (both regular and chaotic) orbits were 
used. 



3. N-BODY MODELS 

We verified that our orbital superpositions correspond to 
true equilibria by realizing them as ./V-body models and in- 
tegrating them forward in time in the gravitational poten- 
tial computed from the N bodies themselves and from the 
point mass representing the black hole. Initial conditions 
were prepared by re-integrating all orbits with nonzero oc- 
cupation numbers. The position and velocity of each tra- 
jectory were recorded at fixed time intervals; the length of 
each integration was determined by the orbit's occupation 
number. The sense of rotation of the tube orbits was cho- 
sen randomly such that the mean motion was everywhere 
zero. The number of particles was 2.42 x 10 6 for 7 = 1 
and 2.19 x 10 5 for 7 = 2; a smaller N was chosen for 7 = 2 
since the ./V-body integrations were slower for the more 
condensed model. 

The initial conditions were then advanced in time us- 
ing th e TV-body code GADGET flSpringel, Yoshida fc White- 



2001 ) , a parallel tree code with variable time steps. The in- 
terparticle softening length was chosen to be 0.005 (0.003) 
for 7 = 1 (2). Figure 2 shows the evolution of the axis 
ratios of the models, computed using the iterative proce- 
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dure described by Dubinski & Carlberg (1991). The axis 
ratios fluctuate with time, but the models retain their tri- 
axial figures very well over many dynamical times. No 
signficant evolution was seen in the radial distribution of 
matter. 




5. DISCUSSION 

We have shown that long-lived triaxial configurations 
are possible for nuclei containing black holes, whether con- 
structed purely from regular orbits, or from a combination 
of regular and chaotic orbits. The latter solutions are re- 
markable given the large fraction of the mass on chaotic 
orbits, of order 50% or more (Table 2). Such solutions 



violate Jeans's theorem in its standard form (e.g. Bin 



ney & Tremaine 1987) but are consi stent with a gener 
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Fig. 3. — Line-of-sight velocities of the rotating regular models 
described in the text, as seen along the y-axis. Contours of positive 
(negative) velocity are represented by full (dotted) lines. 



4. OBSERVABLE PROPERTIES 

Two observational signatures of triaxiality are isophotc 
twists , and kinematic misalignments. Isophotc twists re- 
quirc a variation with radius ot the intrinsic axis ratios, 
which is not a feature of our models. Kinematic misalign- 
ments occur whenever the intrinsic angular momentum 
vector does not coincide with the minor axis of the figure, 
and are generic features of triaxial systems when b oth the 



alized Jeans's theorem (Merritt 1999) if we assume that 
the chaotic building blocks are "fully mixed," that is, that 
they approximate a uniform population of the accessible 
phase space. This appears to be the case based on visual 
inspection of the time-averaged chaotic trajectories. While 
a sudden onset of chaos can effectively destroy triaxiality 
in m odels containing a large p o pulation of regu lar box or- 
bits ( |Merritt fe Quinlan 199§| ; |Sellwood 200l| ), our work 
shows that at least the central parts of galaxies containing 
black holes can remain triaxial even when dominated by 
chaotic orbits. 

Galaxy nuclei are typically modelled as oblate axisym- 
metric systems when derivin g black hole masses from 



stellar-kincmatical data (e.g. van der Marel et al. 1998 



Joseph et al. 2001). Modelling a triaxial nucleus as ax- 
isymmetric will generally lead to biassed estimates of M^h] 
for instance, if an elongated bar is viewed end-on, the ve- 
locity field can mimic that produced by a black hole even 
in th e absence of a central mass concentration ( Gerhard] 
1988). The high-resolution data on which such mass es- 



short- and long- axis tube orbi ts are populated (e.g. Franx 



Illingworth fc de Zeeuw 1991 ) 



timates are based are usually obtained from single slits, 
making it difficult to rule out triaxial shapes using infor- 
mation like that in Figure 3. 

Axisymmetry is typically also assumed when deriving 
rates of interaction o f stars with nuclear black hole s (e.g. 
Syer fc Ulmer 1999|; pVlagorrian fc Tremaine 1999|). The 



We illustrate the possibility of kinematic misalignments 
in our models in Figure 3. Line-of-sight velocity contours 
are plotted for regular models after a fraction of the x- 
and z-tubes were reversed, giving the models net rotation 
about both the long and short axes. For 7 = 1, all x-tubes 
in Figure 3 rotate in the same sense while 60% of the z- 
tubes at each energy rotate in one sense and the remaining 
40% in the other. For 7 = 2, all x-tubes rotate in the same 
sense around the x axis and all z-tubes rotate in the same 
sense around the z-axis. Rotation of the projected models 
along the apparent minor axis of the figure ("minor-axis 
rotation") is apparent. 



triaxial models which we constructed from regular orbits 
were dominated by z-tubes (Table 2), similar to the tube 
orbits that populate an oblate spheroid. The rate of cap- 
ture or destruction of stars by the black hole in these 
models would be similar to the rates in an axisymmet- 
ric nucleus. However the dominance of chaotic orbits in 
our second set of solutions implies a qualitatively different 
sort of behavior: the "loss cone" of stars captured by the 
black hole would be refilled on a time scale determined by 
the frequency of near-center passages of the chaotic orbits, 
rather than by two-body scattering. 

This work was supported by NSF grants AST 96- 
17088 and 00-71099 and by NASA grants NAG5-6037 and 
NAG5-9046 to DM. 
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